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The nonlinear evolution of the m = 1 internal kink mode is studied numerically in a setting where 
the tokamak core plasma is surrounded by a turbulent region with low magnetic shear. As a starting 
point we choose configurations with three nearby g = 1 surfaces where triple tearing modes (TTMs) 
with high poloidal mode numbers m are unstable. While the amplitudes are still small, the fast 
growing high-m TTMs enhance the growth of the m = 1 instability. This is interpreted as a fast 
sawtooth trigger mechanism. The TTMs lead to a partial collapse, leaving behind a turbulent belt 
with g ~ 1 around the unreconnected core plasma. Although, full reconnection can occur if the 
core displacement grows large enough, it is shown that the turbulence may actively prevent further 
reconnection. This is qualitatively similar to experimentally observed partial sawtooth crashes with 
post-cursor oscillations due to a saturated internal kink. 



I. INTRODUCTION 

Abrupt ejection of thermal energy and particles from 
a magnetized high-temperature plasma is frequently ob- 
served in astrophysical and laboratory plasmas. In the 
case of tokamak experiments, internal disruption events 
known as sawteeth bring about a sudden collapse of the 
core temperature through an internal kink instability. A 
thorough understanding of the underlying physical pro- 
cesses is important for the efficient operation of a toka- 
mak device and the prospective application of the toka- 
mak concept for thermonuclear fusion reactors. 

The aim of this paper is to demonstrate the effect of 
magnetohydrodynamic (MHD) turbulence on the evolu- 
tion of the m ~ 1 internal kink mode. Motivated by 
recent results on the instability of current-driven high-m 
multiple tearing modes P, 01 we choose configurations 
with three q = I resonant surfaces located a small dis- 
tance apart. Here, q is the tokamak safety factor (mea- 
suring the field line pitch) and its central value is taken 
to be well below unity {qo < 1). The plasma is taken to 
have a finite resistivity to enable magnetic reconnection. 
This system is, in addition to the resistive m = 1 inter- 
nal kink mode, unstable to a broad spectrum of triple 
tearing modes (TTMs) with helicity q^cs = m/n = 1. At 
present, we neglect two-fluid and kinetic effects and also 
ignore the roles of finite pressure and toroidal curvature. 

In this setting we address some open questions with 
regard to internal disruptions in tokamaks 0, 0, IS 
These include the issue of the rapid onset of a sawtooth 
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crash, known as the trigger problem, and the possibility 
of partial reconnection and compound sawtooth crashes 
(e-g-, QilMI^)' The fast sawtooth trigger is defined as a 
sudden transition from slow growth or stability to rapid 
growth of the to = 1 mode (e.g., 0,0|). The partial 
collapse is defined as a sawtooth crash during which the 
central core region remains intact, so that go remains 
below unity (e.g., 0,0|). It is clear that explanations 
for the sawtooth trigger and partial reconnection events 
require nonlinear effects. Several possible mechanisms 
have been proposed in the past (see, e.g., Ref. for a 
review). Of particular interest here are scenarios that 
consider dynamics related to plasma turbulence. 

In regimes where the to = 1 mode amplitude is still 
small, high-771 modes were previously shown to affect 
the gro wth of the to = 1 instability. Micro-turbulence 
fl5l | (oscillating modes) and TTMs (purely grow- 
ing modes) were found to enhance the m = 1 growth rate. 
In other related studies it was shown that the to = 1 
mode can be stabilized by local oscillations at the reso- 
nant surface 0, 0| and micro-turbulence in the region 
of the thin current layer may lead to enhanced effective 
resistivity and thus faster reconnection 0, ITgl l20l l2ll| . 
Conversely, viscosity may be increased, which reduces 
the reconnection rate |23|. In a broader sense, mag- 
netic braiding and field line stochasticity may also be 
viewed as "turbulent" dynamics, and these were also 
shown to yield enhanced growth rates for reconnecting 
modes [2J, |2J| . With regard to the long-time nonlinear 
evolution, pressure-driven instabilities (e.g., ballooning) 
were shown to be able to lead to a saturation of the inter- 
nal kink at finite amplitude (partial collapse) p5ll2^l27l |. 

The scenario considered here is comparatively simple 
and includes only a minimum of physical effects. In the 
first part of thispaper the case of a TTM-driven inter- 
nal kink mode [l| is examined and two new results are 
presented: establishment of a fast growing to = 1 mode 
structure at low amplitudes and, in regimes with high vis- 
cosity {Pr ^10), explosive growth during the transition 
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to the turbulent regime. 

The transition to turbulence occurs via a partial (an- 
nular) collapse, whereupon a turbulent belt forms around 
the central core. Indeed, in Ref. the conjecture was 
made that partial sawtooth crashes may be associated 
with wide-spread MHD turbulence in the reconnected re- 
gion. In Ref. Q the possibility for the m = 1 mode to 
saturate in such a state was demonstrated using a nu- 
merical simulation where the core was constrained to a 
linear motion in the poloidal plane. 

The second part of this paper deals with the long-term 
evolution of the internal kink, while surrounded by a low- 
shear region with g ~ 1, governed by MHD turbulence. 
The scenario considered here is more generic than that in 
Ref. since we apply a fully random perturbation and 
allow each Fourier mode to alter its poloidal phase an- 
gle through interaction with other modes. The resulting 
changes in the kink flow lead to an irregular or "mean- 
dering" motion of the core. One case is presented where 
full rcconncction occurs eventually. In another case, the 
kink is found to saturate and decay, which shows that 
the MHD turbulence may prevent full reconnection. 

A reduced model is used for the simulations in order to 
obtain first insights on a fundamental level and to lay the 
foundations for further investigations with more realistic 
models. Physical effects ignored in this work, such as 
two-fluid, curvature, and finite-beta effects, may play a 
significant role. In future work it would be interesting to 
see under which conditions and in which way these effects 
alter certain quantitative and qualitative features of the 
results presented here, including the linear and nonlinear 
instability growth, the dominant mode numbers in the 
inter-resonance region, and the evolution of the internal 
kink. 

This paper is organized as follows. The physical model 
is introduced in Sec. HTIand the numerical method is de- 
scribed in Sec, mil The equilibrium used, its linear insta- 
bility characteristics and the initial perturbation applied 
are given in Sec. IIVI In Sec. |V] we present results on 
the early evolution of the m = I mode in the presence 
of fast growing high-m TTMs. The transition to turbu- 
lence is treated in Sec. IVII and the long-term evolution 
is described in Sec. IVIII A summary, further discussions 
and conclusions are given in Section IVIIII 



II. MODEL 

We use the reduced set of magnetohydrodynamic 
(RMHD) equations in a cylindrical geometry in the limit 
of zero beta [2^ Is^l . This model has proven to be use- 
ful in studies of MHD instabilities when the focus is on 
a qualitative description of fundamental aspects of the 
magnetized plasma system, as is the case here. The 
RMHD model governs the evolution of the magnetic flux 
function ip and the electrostatic potential (/), as described 



in Ref. . The normalized RMHD equations are 

M = 0] - ac<^ - ^Hp iVJ ~Eo) (I) 
dtu = [u, (j)] + [j, i'] + da + Re^l^lu. (2) 

The time is measured in units of the poloidal Alfven time 
''Hp = ^/JkiPnia/ Bq and the radial coordinate is normal- 
ized by the minor radius a of the plasma. Here, p,,! is 
the mass density and Bq the strong axial magnetic field. 
The current density j and the vorticity u are related to 
ip and (j> through j = — V^V ^-^d = respectively. 

In order to provide a simple mechanism for mag- 
netic rcconncction, a resistive diffusion term is included 
in Eq. Its strength is measured by the magnetic 

Reynolds number ^Hp = Tij/thp, with = a^fio/rjo be- 
ing the resistive diffusion time and rjo = ri{r ~ 0) the elec- 
trical resistivity in the plasma core. We use ^Hp = 10^, 
which is numerically efficient and physically reasonable 
in the framework of the model used. Flow damping is 
provided by an ion viscosity term in Eq. Q . Viscous dis- 
sipation is measured by the kinematic Reynolds number 
Reup = a^/i/THp, where v is the ion viscosity. Long-time 
calculations are performed for i?eHp — 10^ and 10^, as 
will be specified case by case. 

In order to ensure that, in the absence of magnetic 
reconnection, the equilibrium remains unchanged, the 
source term S^^Eq is included in Eq. (Q. With Eq = fjj 
it compensates the resistive dissipation of the equilib- 
rium current. The loop voltage measured by i?o is 
taken to be constant, so the resistivity profile is given 
in terms of the equilibrium current density distribution 
as f}{r) = j(r = 0)/j(r). For simplicity, the temporal 
variation of the resistivity profile fj is neglected. 

As in Ref. each field variable / is decomposed into 
an equilibrium part / and a perturbation / as 

/(r,,?,C,t)=7(r) + 7(r,^,C,i). (3) 

The system is described in terms of the Fourier modes, 
')prn.,n and (j)m,m obtained from the expansion 

fir, §,Ct)^lj2 fr-Ar, t).e^(™''-"'^) + c.c, (4) 

rn.n 

with m being the poloidal mode number and n the 
toroidal mode number. In the following, the (m, n) sub- 
scripts will often be omitted for convenience. We consider 
only the dynamics within a given helicity h ^ m/n = 
const, so the problem is reduced to two dimensions. 

For the description of the dynamics in this system, 
it is useful to define the helical flux function -0* with 
corresponding current density as 

71 2?^ 

ip*^'ip + —r'^ and = = j . (5) 

Zm m 

The evolution of the individual Fourier modes is de- 
scribed in terms of their kinetic and magnetic energies, 

i^™"„ = |V0™,„P and S-^'s^lVV'™,™^ (6) 
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and the corresponding nonlinear growth rates, 
dln£;i=;" 



/m,n 



it) 



and YZSit) 



d In E""^s 



(7) 



2dt 2di 

(these are ampHtude growth rates, hence the factor 1/2) 
1 

In Eq. ©, |/m,«P = / dr.r.Cm|/m,„(r)|2, with Cm=o = 



47r and Cm^o ~ 27r. 



III. NUMERICAL METHOD 

For the numerical solution of the model equations 
and a two-step predictor-corrector method is applied. 
In the first time step the dissipation terms are treated 
implicitly, all others explicitly, and the field variables are 
estimated at an intermediate time step t + At/2. The 
second is a full time step, t t + At, with the right-hand 
sides of Eqs. ^ and ^ evaluated at the intermediate 
time step t + At/ 2 estimated before. In the nonlinear 
regime the time step size is of the order At ~ 10^^. 

128 Fourier modes (including to = 0) are carried. 
The Poisson brackets [f,g] = ■^{drfd^g — drgd^f) are 
evaluated in real space. This pseudo-spectral method 
has been applied together with an appropriate dealiasing 
technique. The outcomes of the long-term evolution in 
both cases studied, full reconnection and kink saturation, 
were also confirmed using 256 Fourier modes. 

The radial coordinate is discretized using a non- 
uniformly spaced grid, with a grid density of up to 
A^,r^ « 1/2000 [1/5000 for numerical checks] in regions 
where sharp current density peaks occur. A fourth-order 
centered-finitc-diffcrence method is applied for the dr- 
tcrms in the Poisson brackets. The Laplacians , > ~ 

^drrdr — rr? jr^ are evaluated at second-order accuracy 
(tridiagonal matrix equations). 

Periodic boundary conditions are applied in the az- 
imuthal and axial directions. At r = 1 an ideally con- 
ducting wall is assumed, requiring all perturbations to 
be identical to zero at that location: f{r = 1) = (fixed 
boundary, no vacuum region). At r = 0, extraneous 
boundary conditions are applied to ensure smoothness: 
drI„i=o{r = 0) = and fm^r = 0) = 0. 



IV. EQUILIBRIUM, LINEAR INSTABILITY 
AND INITIAL PERTURBATION 
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FIG. 1: (Color online). Safety factor profiles q{r). Resonant 
surfaces are indicated by circles. Model parameters for both 
cases, labeled (T-1) and (T-2), are given in Table 13 The 
profile properties are listed in Table [H] and the dispersion 
relations are shown in Figs. Inland |21 
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FIG. 2: (Color online). Spectra 7iin(rn,) of unstable eigen- 
modes in Case (T-1) for Shp = 10^ and ReHp = lO''. The 
growth rates 7ii„ , 7ii^' and yli^ of the three eigenmodes M*-^-*, 
M(2) and M(3) (cf. Fig. 4 in Ref. Q) are shown. 



flux function and current density profiles are obtained 
though the relations 



q ' = ---j-i^o,o and jo,o 
r ar 

In this paper configurations with three resonant sur- 



1 d 

. (9) 

r ar q 



The equilibrium state is taken to be axisymmetric 
(only TO = n = components) and free of flows, i.e., 

= u = O. (8) 

The equilibrium magnetic configuration is uniquely de- 
fined in terms of the safety factor q(r), and the magnetic 



Case 


go 


TA 


MO 


Ml 


m 


n 


h 


ru 


ri2 


(T-1) 


0.73 


0.455 


0.93 


1.45 


1 


1 


-0.09 


0.5406 


0.039 


(T-2) 


0.73 


0.455 


0.93 


1.45 


1 


1 


-0.098 


0.5666 


0.064 



TABLE I: Parameter values for the q profiles shown in Fig.0 
using model formula (11) in Ref. 0. 
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FIG. 3: (Color online). Spectra ^unim) of unstable eigen- 
modes in Case (T-2) for Shp = 10® and i?eHp = 10*. The 



(1) 



"(lin '^^ three eigenmodes M'^-*, 



growth rates 7]j 

Af(^> and Af^^> (of. Fig. 4 in Ref. 0) are shown. 



Case 


Di2 


D23 


Sl 


S2 
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Vi- 
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V3 




(T-1) 


0.053 


0.043 


0.35 


-0.56 


1.20 


1.7 


1.1 


3.5 


13 


(T-2) 


0.052 


0.052 


0.23 


-0.23 


0.59 


1.6 


1.2 


1.9 


13 



TABLE II: Properties of the q profiles shown in Fig. Q The 
mode number of the fastest growing mode, r/ipeak (cf. Figs.|21 
andOJ, is valid for Sup = lO*' and J?eHp = 10® [Case (T-1)], 



iJcHp = 10** [Case (T-2)]. 



faces with qs = q{rsi) = m/n ~ 1, which arc located 
at radii Vsi {i = 1,2,3), are considered. The distances 
between the resonances, Dij = \rsj —rsi\, are chosen suf- 
ficiently small, so that the spectra of unstable TTMs are 
broad and the dominant modes have m ^ 0(10). The 
equilibrium q{r) profiles used arc shown in Fig. ^ They 
are obtained using the model formula given by Eq. (11) 
in Ref. 01 with the parameters in Table The two cases 
studied are labeled (T-1) and (T-2). 

The dispersion relations (spectra of linear growth 
rates) 7iin("T-) for all unstable eigenmodes are given in 
Fig. 121 for Case (T-1) and in Fig. 01 for Case (T-2). The 
linear eigenmode structures for Case (T-1) were shown 
in Ref. _.l,] and arc similar for Case (T-2). Let us recall 

that eigenmode M^^^ (with growth rate jI^^}) is associ- 
ated with the resonant surface r = rsi and is an ordinary 
(single) m = 1 internal kink-tearing mode (stable for 

771 > 1). A/^^) (with 7!;^^) is associated with in the 
sense that it is active in the region < r < for 777 = 1 
and Tgi < 7- < for 777 > 1. Thus, it may be regarded 
as a double tearing mode (DTM). Finally, M^^) (with 

/ON 

Tiin ) associated with and may be regarded as the 
actual TTM eigenmode. In Figs. Inland O the dominant 
eigenmodes and the 777 = 1 mode arc indicated by arrows 
labeled with the corresponding mode numbers. 

The characteristics of the equilibrium configurations 
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FIG. 4: (Color online). Early evolution in Case (T-1). (a): 
Kinetic energies [Eq. Q] of the m — 1 mode and the 

two fastest growing modes m = 13 and m — 14. (b) and (c): 
Magnetic and kinetic growth rates 'j^'^n and 7m"n [Eq. Q]. 
The main stages to be distinguished are: (i) linear growth, (ii) 
nonlinearly driven growth, (iii) transition to the fully nonlin- 
ear (turbulent) regime. The value of the driven growth rate 
expected during stage (ii) is 7drivc ~ 27iin(mpoak) ~ 16x10"^, 
as indicated in (b). Sup = 10®, i?eHp = 10®. 



are summarized in Tabled including the values of the 
resistivity profile r) at the resonances and the poloidal 
mode number 777poak of the dominant TTM eigenmode 
obtained from Figs. |2 and |2| We will see that the value 
of 777,poak is an important clue for the interpretation of the 
nonlinear dynamics since it dictates a structure size that 
is most likely to be encountered in the poloidal direction. 

Starting from an unstable equilibrium, the instability 
is excited by applying an initial perturbation of the form 



1 

^(^ = 0) = 9 E *o,™r(r - l)e 



(10) 



where ^'o.m = 10^ is the perturbation amplitude and 
I?* = 7? — q^^C is a helical angle coordinate. Each mode 
has an initial poloidal phase shift 7?o,m with a randomly 
assigned value in the range 79o.m G [0, tt]. 
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FIG. 5: (Color online). EflFect of the nonlinear driving on the 
m = 1 mode structure. In (a) the linear mode structure of the 
fastest growing m = 1 eigenmode, V'l'hi ("^ = l)i is shown. The 



two fastest growing modes, nipcak ~ 13 and mpcak + 1 = 14 
(cf. Fig. 121 , give rise to a rapidly growing radially localized 
m — 1 component through the convective nonlinearity 
in Eq. Q . In (b) and (d) the structures of the two main terms 
are plotted [(b) dominates]. In (d) we show the situation 
before (t = 120) and after the new fast growing m — 1 mode 
structure is established {t — 140 and 180). The curves in (a) 
and (d) are scaled such that the amplitude of the internal kink 
component (0 < r < Vsi) is the same at all times. Resonant 
surfaces are indicated by dotted vertical lines. 



V. EARLY DYNAMICS: ESTABLISHMENT OF 
NONLINEAR MODE STRUCTURE AND FAST 
GROWTH OF THE m = 1 MODE 



nantly mpoak = 13 and mpeak ± 1- Since 7iin(mpeak) ~ 
7iin("^pcak±l) the driven growth rate of the m = 1 mode 
is 7drivc ~ 27iin(mpcak) « 16 X 10"^. The purely nonlin- 
early driven m = mode grows at the same rate, as can 
be seen in Fig.^fb) [phases (i) and (ii)]. Note that in the 
present case 7drivc this is almost one order of magnitude 
higher than the linear growth rate of the m = 1 mode, 
7iin(m = 1) = 1.7 X 10-2. 

The effect of the nonlinear driving on the m = 1 mode 
structure can be seen in Fig. |31 In Fig. E^a) the hnear 
mode structure of the TTM-type eigenmode is plotted, 

specifically the flux function ipuni''^ ~ "^^^ shapes 
of the driving terms arising from the convective nonlin- 
earity [V', in Eq. |^ arc plotted in Fig. [Sfb) and (c). 
In Fig. Eld) it can be seen how this radially localized 
driving appears in the nonlinear m = 1 mode structure 
'0drivc(™ = 1) in the form of "spikes" located near the res- 
onant surfaces. It is particularly remarkable that, after a 
certain ratio between the amplitudes of the driving com- 
ponent (around Tss) and the component remaining from 
the linear eigenmode structure (in the region < r < Vgi) 
is reached, the mode structure docs not change further 
[Fig. Eld), t = 140-160] and grows as a whole at the en- 
hanced rate 7drive- Thus, a new fast growing global m = 1 
mode is created through radially localized driving. Com- 
parisons made with other simulation runs indicate that 
the relative size of the driving component and the peak 
in the region < r < Tgi in the established nonlinear 
mode structure depends on the ratio between 7drivc and 
7iin(m = 1): If 7drivo/7iin(m 1) IS large, the driving 
component (around rg^) has a large amplitude compared 
to the kink component (0 < r < rgi), and vice versa. 



Due to the disparate growth rates between the m = 1 
mode and the fastest growing TTMs [here, 7peak ^ 
5 X 7ii„(m = 1); cf. Figs. |21 and |3] nonlinear interac- 
tions begin already in regimes where mode amplitudes 
are still small. In this section we focus on these early 
stages of evolution where turbulence has not yet devel- 
oped. They are most conveniently studied by consider- 
ing the dynamics of individual Fourier modes. Results 
are presented only for Case (T-1), using Sup = 10^ and 
^eHp = 10^. Similar behavior is observed in Case (T-2). 

Figure ^a) shows the evolution of the kinetic energies 
of the TO = 1 mode and the two fastest growing modes 
TOpeak = 13 and TO = 14 (cf. Fig. The corresponding 
growth rates are shown in Fig. 01b) and (c) . It can be 
seen that the to = 13 and to = 14 modes grow at their 
linear growth rates until t w 170 and saturate during 
170 ^ t < 200. The to = 1 mode grows linearly until 
t w 100 [phase (i) in Fig. 21- Subsequently, its growth 
rate increases and exponential growth continues at an 
enhanced rate during 100 < t < 170 [phase (ii)]. Upon 
entering the nonlinear regime, the to = 1 growth rate 
drops [phase (iii)]. 

The enhanced growth during phase (ii) is due to non- 
linear driving by the fastest growing modes, predomi- 



VI. BEHAVIOR OF THE m = 1 MODE DURING 
THE TRANSITION TO TURBULENCE 

Let us now consider the transition to the fully non- 
linear, turbulent regime. In Fig. 0] this transition takes 
place via a gradual decrease in the growth rates dur- 
ing phase (iii). However, it turns out that this is only 
the case in regimes where the effect of viscosity is neg- 
ligible. In Fig. El the evolution of the nonlinear growth 
rates 7™i^, 7i.'i and 7™o^ is shown for magnetic Reynolds 
numbers 5hp = 10^, 10^, 10*^ (from left to right) and the 
kinematic Reynolds numbers i?eHp = 10^, 10^, 10^ (from 
top to bottom). Thus, the Prandtl number Pr, which 
measures the relative strengths of viscosity and (resis- 
tive) diffusion (Pr = Sup/Reup oc v/rjo) varies over the 
range 10" ^ < Pr < 10^ (bottom left to top right). The 
growth rate of the to = mode, 7™o^, is shown since it 
clearly reflects the level of nonlinear driving at all times. 
Qualitatively, the results in Fig. El niay be summarized 
as follows: 

1. Pr ^ 0.1 [Fig. Elg)]: During the transition to the 
nonlinear regime, 7^^™ first decreases slowly {t sa 
150-180) and then drops rapidly. Similar behavior 
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FIG. 6: (Color online). Effect of varying Sup and -Renp on the evolution of the m = 1 mode, described in terms of 71 'i (dashed 
line) and 7™i^ (solid line). The growth rate of the m — mode is plotted as well (dash-dotted line), showing the level of 
nonlinear driving starting at t = 0. We distinguish between linear growth (7iin) and nonlinearly driven growth due to the 
fastest growing modes (7drivc). For Pr > 10 explosive growth follows, up to a maximum growth rate 7max, as indicated in (b). 
The time interval shown is from t = 50 (after the transient relaxation) until the annular collapse and onset of MHD turbulence. 
The configuration used is Case (T-1). The data for Pr = 1 in (d) is the same as in Fig. Qlb) above and in Fig. 5 of Ref. 0. 



is observed for Pr = 10 ^ and is also expected for 
Pr < 10^2 

2. Pr ^ 1 [Fig. |S^d,h)]: Compared to the results ob- 
tained with lower Pr, 7drivc is now reduced due to 
the stabilizingeffect of the viscosity on 7iin("^) (cf. 
Fig. 5 in Ref. During the transition to the fully 
nonlinear regime 7^ remains high for a certain pe- 
riod of time and may exhibit oscillatory behavior 
as in Fig. Efh) before dropping. 

3. Pr - 10,10^,103 [Fig. El;a,b,c,c,f,i)]: The growth 
rate 7drivo is further decreased due to higher vis- 
cosity. However, at the end of the driving phase a 
significant increase in the growth rate from 7drive ~ 
27pcak to a value 7max [indicated in Fig. Elb)] is 
observed. This effect is most pronounced in the 
kinetic growth rate 7^'". Note in particular that 
7max can get close to the value of 7drive obtained 
for Pr ~ 1 [compare, e.g., Fig. |Bfb), (e) and (h)]. 
It is likely that this behavior can also be observed 
for Pr > 10^. 

Let us note that during this rapid growth the nonlinear 
interactions already include many Fourier modes. We 
suspect that the explosive growth phase observed here for 
Pr > 10 is due to these nonlinear interactions becoming 



more important than viscous damping. The latter had 
reduced 7drivc through a reduction of the linear growth 
rates juni'm)- 



VII. ANNULAR COLLAPSE AND EFFECT OF 
MHD TURBULENCE ON THE INTERNAL KINK 

In this section we investigate the long-term evolution 
in Cases (T-1) and (T-2) (Fig.[TJ. We describe the mag- 
netic and E X B flow structures generated through TTM 
reconnection and analyze the evolution of the m = 1 in- 
ternal kink mode while it is surrounded by MHD turbu- 
lence. The values for the dissipation parameters in Case 
(T-1) are 5hp = 10^ and i?eHp = lO'^, as in Sec.|V| For 
Case (T-2) we choose Shp = 10^ and i?eHp = 10^ 



A. Case (T-1): Annular collapse and full 
reconnection 

We begin with a discussion of snapshots taken in Case 
(T-1), which are shown in Figs. [7| and |S1 and labeled (A)- 
(F). The initial perturbation is sufficiently random so 
that reconnection occurs all around the core [Fig. [Tf A)] . 
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FIG. 7: (Color online). Annular collapse and long-term evolution in Case (T-1). The first three snapshots taken at (A) t — 200, 
(B) t = 240, and (C) t = 280 are shown (continued in Fig. Each snapshot consists of contour plots of the helical flux 
(left) and the electrostatic potential cj) (right). Arrows indicate the flow directions. The dashed circles indicate the outermost 
ijj contour of the core and have been superimposed on the contours for clarity. The small diagrams in the middle show the 
instantaneous profiles q{r,t) and (j,) = [j»(r', t)]o,o. -^Hp ~ 10^, Reap ~ 10^. 



However, it is not isotropic, so that some magnetic is- 
lands grow faster than others. The sizes of the islands 
reflect the shape of the spectrum of linear growth rates 
"fiini'm) with dominant modes having m ~ mpoak = 13 
[cf. Fig. 12] . During this annular collapse the amplitude 
of the m = 1 mode is still small and the q{r) profile 



is flattened annularly [Fig. [7fA)]. The magnetic islands 
in the inter-resonance region exhibit complicated coales- 
cence dynamics and a turbulent annular region is created 
around the core [Fig. EfB) and (C)]. In the meantime, 
the core displacement becomes observable. The core tra- 
verses the turbulent belt [Fig. ^D) and (E)] and upon 
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making contact with the flux surfaces beyond the outer- 
most resonant radius (r = rss) core reconnection begins 
[Fig. IHE) and (F)]. In the present case, (T-1), core re- 
connection is hkely to proceed to completion, i.e., full 
reconnection is expected. For completeness, the evolu- 
tion of the m — I energies and growth rates is shown in 
Fig. Eta) and (b), respectively. 

An important observation is the following. While the 
TO = 1 mode was originally perturbed with i9o.m=i = 



(i.e., core motion in the positive x direction), the kink 
flow continuously changes its direction, as is obvious from 
the arrows drawn along with the (j) contours in Figs. [7| 
andlHI Note that the core does not rotate; merely the 
direction of its translational motion alters (in RMHD, 
(/>o,o = at all times if it is zero initially). The core's 
motion is quantified and shown in more detafl in Fig.El^c), 
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FIG. 9: (Color online). Long-term evolution of the m = 1 
mode in Case (T-1). (a); Magnetic and kinetic energies E^'!^ 
and . (b): Magnetic and kinetic growth rates 7™i^ and 
7i,'i • (c): Motion of the core in the kink flow with dominant 
mode number m = 1. The magnetic axis is located at the 
poloidal angle iJmag and is moving in the direction i^kSU'' (in 
the poloidal plane at = 0) [see deflnitions in Eq. HlU l. In 
(a), (b) and (c) the times at which snapshots were taken are 
indicated by circles and labeled (A)-(F) [cf. Figs. Q and 



using the angles t^J^Iag s-'^d i^^Jn^^- These are defined as 
J dr.rIm{V'i.i(?')} 



/ d?'.rRc{V'ia(?')} 
/ dr.rRe{0i,i(r)} 
J dr.rlni{(/ii,i (r)} 



(11a) 
(lib) 



(integration interval: < r < 0.07) and give an approx- 
imate image of the core's motion when the kink ampli- 
tude is not too large. During the linear phase (i) both 
angles arc zero, in agreement with the initial perturba- 
tion. During the driving phase (ii) the angles switch to 
^mag ~ '^kin'^ ~ 3° (determined by the driving modes). 
After the transition to turbulence in phase (iii) . the direc- 
tion of the kink flow varies frequently, so that the core's 
motion becomes rather complicated. 



B. Case (T-2): Kink saturation and partial 
reconnection 

Here we discuss the long-term evolution in Case (T-2). 
Snapshots are presented in Figs.llOlandllll In Fig.ll2lthe 
evolution of (a) the m = 1 energies, (b) corresponding 
growth rates, and (c) the core's location and kink flow 
direction are shown. Note that the early evolution is 
very similar to that of Case (T-1) described in Section 
IVl It consists of (i) linear growth, (ii) nonlincarly driven 
fast growth, and (iii) transition to the turbulent regime 
with gradually decreasing growth rates (here Pr = 0.01). 
We may thus omit the further discussion of these stages, 
referring to Sec. Ivl above. 

Similarly to Case (T-1), the first macroscopically ob- 
servable event is an annular collapse due to high-m TTMs 
without significant displacement of the core [Fig. lTUT A)]. 
Again, the q{r) profile is flattened in the inter-resonance 
region. Subsequently the m = 1 mode grows inside the 
turbulent belt [Fig. llOf B)]. However, in contrast to Case 
(T-1), here the m — 1 mode saturates after reaching a 
relatively large amplitude [Fig. lTUT C)]. This occurs att ~ 
962, as can clearly be seen in the m = 1 magnetic energy, 
and the associated growth rate, 7™f® [Fig. IT^ a) 
and (b)]. The direction of the kink flow reverses, as is ob- 
vious from Snapshots (C)-(E) (Figs.lTUland lllll and from 
the 180-degree jump in "d^l^ in Fig. I12r c). Afterwards, 
the kink amplitude decays [Fig. Ilir E)]. overshoots, and 
grows again in a different direction [Fig. Illf F)]. 

The saturation of the m = 1 mode observed here seems 
to be due to an island-like structure developing "in front" 
of the displaced core, when the latter approaches the pe- 
riphery. The island remains trapped, i.e., it is not ex- 
pelled in the poloidal direction. Consequently core recon- 
nection takes place at two separate points which in the 
present case are located an angle Ad sa 116° apart, as in- 
dicated in Fig. lTUT C). Since each reconnected flux surface 
adds to the island's width, this structure can counter the 
internal kink, induce a rebound and send the core back 
into the center. This scenario is realized in the present 
case. 



C. Modulation of the kink flow 

In addition to the kink flow changing its direction, at 
certain times we observe an m = 1 modulation of the 
contours in the core's interior. In Case (T-1) this occurs 
around t ~ 400, Snapshot (D) [Fig. This modulation 
of the E X B drift velocity is not strong enough to visi- 
bly alter the -0 contours in the core, which therefore re- 
main circular (not shown) . Further details can be seen in 
Fig. El where the profile of the radial velocity Vr oc (p/r 
is shown at several times in the interval 360 < t < 510. 
Some peaks in Vr, like those near the magnetic axis, per- 
form oscillations in the manner of a standing wave [left- 
hand side in Fig. I13f b)]. Other peaks, like those in the 
region 0.2 < r < 0.35, do not change their signs until 
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FIG. 10: (Color online). Annular collapse and long-term evolution in Case (T-2). The first three snapshots taken at (A) 
t = 500, (B) t = 750, and (C) t = 950 are shown (continued in Fig.ITTJ- •S'hp = 10^ Re^p = 10*. 



about t « 500 [Fig. C^c)]. The radial wavelength of VIII. DISCUSSION AND CONCLUSIONS 

the modulation is observed to change relatively slowly. 
Typically, it measures between 1/2 and 1/5 of the core's 

radius. Realizations of this m = 1 modulation can also We have studied the nonlinear evolution of the rn = 1 
be observed in Case (T-2) (Figs. [Hand [ill). internal kink mode in a configuration with three q = 1 

resonant surfaces where high-m TTMs are strongly un- 
stable and lead to an annular collapse. The latter leaves 
behind a turbulent belt around the unrcconnccted core 
plasma. The simulation results show that high-m TTMs 
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and MHD turbulence, which are locahzed in an annular 
region, are able to strongly affect the evolution of the 
TO = 1 internal kink mode, which is a global instability. 
We conclude that multiple tearing modes (here, TTMs) 
and MHD turbulence may play a significant role during 
partial, compound or full sawtooth crashes in tokamak 
plasmas, as will be discussed in the following. 

In the beginning, a fast sawtooth trigger, defined as a 
sudden transition from slow to rapid growth, was real- 



ized: after a phase of slow linear growth, rapidly grow- 
ing q = 1 TTMs give rise to a new fast growing non- 
linear m = 1 mode. This instability reaches an observ- 
able amplitude within a time much shorter than expected 
from the linear growth rate. Moreover, the transition to 
the fully nonlinear (turbulent) regime occurs via a phase 
of explosive growth when the Prandtl number is large 
{Pr > 10). As proposed in Rcf. 1] this enhancement 
of the m = 1 mode due to high-m TTMs is a possible 
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FIG. 12: (Color online). Long-term evolution of the m — 1 
mode in Case (T-2). Arranged as Fig. |5] In (a), (b) and 
(c) the times at which snapshots were taken are indicated by 
circles and labeled (A)-(F) [cf. Figs. ITHlandin]. 

mechanism for the fast sawtooth trigger, provided that 
multiple q = 1 resonant surfaces are formed during the 
sawtooth ramp. 

During the further evolution, the turbulence in the col- 
lapsed annular region was seen to alter the direction of 
the kink flow responsible for the core displacement. Both 
continuing growth [Case (T-1)] and saturation of the kink 
[Case (T-2)] were observed, which shows that full as well 
as partial reconnection may occur in this setting. It was 
also found that the effect of the turbulence was not lim- 
ited to the collapsed annular region and the overall mo- 
tion of the core inside this turbulent belt: Perturbations 
of the electrostatic potential were even found to pene- 
trate into the central core region in the form of an m = 1 
modulation on top of the kink flow. 

The converse effect, i.e., the influence of the core dis- 
placement on the surrounding turbulence has not been 
addressed and is left for future study. This is expected 
to be important since (a) the core displacement changes 
the geometry of the turbulent region and (b) the return 
flows of the internal kink arc likely to interact with the 
turbulence. One particular question to be addressed is 



FIG. 13: (Color online). Evolution of the m = 1 modulation 
on the radial displacement velocity profile ii^ oc 0/r in Case 
(T-1). (a): Complete Vr profile along the poloidal angle = 
15°, i.e., roughly where the magnetic axis is located. The time 
is t — 400 and coincides with Snapshot (D) in Fig. |H| (b): 
Temporal evolution of the m = 1 modulation in the interval 
360 < t < 510 (from top to bottom). The radial location 
of the magnetic axis varies as indicated by the circles [e.g., 
r^Mt = 360) = 1.6 X 10-^ and r^^is{t = 510) = 2.7 x 10""]. 
The average radial velocity of the core is about Vr ~ 0.7 x 10^^ 
(in a/rHp). (c): Same data as in (b) but in the radial interval 
0.1 < r < 0.5 and scaled up for clarity. 

whether and how the core contributes to the formation 
of the trapped island that prevents further reconnection 
in Case (T-2) [cf. Fig.CUC)]. 

Our results agree with some aspects of the partial 
sawtooth crash scenario suggested in earlier studies: A 
"shoulder" on the q{r) profile forms where g ~ 1, and this 
region i s g overned by electromagnetic turbulence (e.g., 
|28L I3lll33 |). Indeed, the conjectures made by these au- 
thors imply that continued growth of the m — \ mode [as 
in Case (T-1)] must be prevented by some means, so that 
the partial collapse remains partial. We have demon- 
strated that MHD turbulence is one possible mechanism 
leading to a saturation of the internal kink [as in Case (T- 
2)]. The residual core displacement may then account for 
the post-cursor oscillation observed experimentally after 
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partial reconnection events (e.g., [ssIIb^ ). 

The fact that the long-term calculations in the two 
cases considered here have different outcomes, namely 
full reconnection in Case (T-1) and partial reconnection 
in Case (T-2), requires commenting on. Note that the 
initial conditions do not differ largely, except for the lin- 
ear growth rates in Case (T-1) being twice as high as in 
Case (T-2) due to a higher magnetic shear (cf. Figs. [2| 
andO. While we were able to demonstrate that MHD 
turbulence provides prerequisites for a saturation of the 
internal kink, the ultimate goal of determining quanti- 
tative criteria for partial reconnection requires further 
investigations using more realistic physical models. As 
mentioned in the introduction, potentially important ef- 
fects to be considered include two-fluid, curvature, and 
finite-beta effects. In particular, if some of these would 
be found to have a stabilizing influence on the high-m 
modes, which are essential in the present work, the re- 
sults may be altered. To our knowledge, such studies 
have not been conducted for a comparable scenario, i.e., 
strongly coupled multiple tearing modes such as DTMs 
or TTMs. 

In this study. MHD turbulence was generated through 
current-driven resistive TTMs, requiring multiple q — 1 
resonant surfaces. The inclusion of a coUisionless recon- 
nection mechanism is expected to yield higher kinetic 
energies jssf and thus stronger turbulence interacting 
with the internal kink. Furthermore, we conjecture that 
our principal result, namely that the internal kink mode 
can be strongly affected by tcaring-modc-driven MHD 
turbulence, will also apply when the latter is generated 
by some other means, such as pressure-driven MHD or 



micro- instabilities [36|. 

Let us note that there are several other mecha- 
nisms which were proposed as possible explanations for 
the rapid collapse and partial crash phenomena, which 
were discovered assuming different pre-crash conditions 
(e.g., Ref. 13 and references therein). However, the cur- 
rently available experimental data is not yet conclusive 
enough to rule out one proposed sawtooth crash scenario 
or another. Moreover, the detailed evolution may vary 
between different machines, shots and even between saw- 
tooth crashes of a single discharge. 

Through recent progress in plasma diagnostics it has 
become possible to detect high-m magnetic islands as- 
sociated with low-order q = m/n resonant surfaces j37l |. 
Thus, experimental checks of our simulation results seem 
feasible in the near future. 



Acknowledgments 

A.B. would like to thank Y. Kishimoto, Y. Naka- 
mura and F. Sano for valuable discussions. He also ac- 
knowledges the Max-Planck-Institut fiir Plasmaphysik in 
Garching for its hospitality and for providing computa- 
tional resources for numerical checks. S.B. acknowledges 
the Graduate School of Energy Science at Kyoto Univer- 
sity and the Center for Atomic and Molecular Technolo- 
gies at Osaka University for their support and hospitality. 
This work was partially supported by the 21st Century 
COE Program at Kyoto University. 



[1] A. Bierwage, S. Hamaguchi, M. Wakatani, S. Benkadda, 
and X. Leoncini, Phys. Rev. Lett. 94, 065001 (2005). 

[2] A. Bierwage, S. Benkadda, S. Hamaguchi, and 
M. Wakatani, Phys. Plasmas 12, 082504 (2005). 

[3] B. N. Kuvshinov and P. V. Savrukhin, Sov. J. Plasma 
Phys. 16, 353 (1990). 

[4] S. Migliuolo, Nucl. Fusion 33, 1721 (1993). 

[5] R. J. Hastie, Astrophys. Space Sci. 256, 177 (1998). 

[6] F. Porcelli, S. Annibaldi, D. Borgogno, P. Buratti, F. Cal- 
ifano, R. Coelho, E. Giovannozzi, D. Grasso, E. Lazzaro, 
F. Pegoraro, et al., Nucl. Fusion 44, 362 (2004). 

[7] H. Soltwisch, W. Stodiek, J. Manickam, and J. Schliiter, 
in Plasma Physics and Controlled Nuclear Fusion Re- 
search 1986 (International Atomic Energy Agency, Vi- 
enna, 1987), vol. 1, p. 263, IAEA-CN-47/A-V-1. 

[8] A. W. Edwards, D. J. Campbell, W. W. Engelhardt, H.- 
U. Fahrbach, R. D. Gill, R. S. Granetz, S. Tsuji, B. J. D. 
Tubbing, A. Weller, J. Wesson, et al., Phys. Rev. Lett. 
57, 210 (1986). 

[9] F. M. Levinton, S. H. Batha, S. H. Yamada, and M. C. 
ZarnstorfT, Phys. Fluids B 5, 2554 (1993). 
[10] A. Y. Aydemir, Phys. Fluids B 4, 3469 (1992). 
[11] X. Wang and A. Bhattacharjee, Phys. Rev. Lett. 70, 1627 
(1993). 



[12] X. Wang and A. Bhattacharjee, Phys. Plasmas 2, 171 
(1995). 

[13] D. Biskamp and T. Sato, Phys. Plasmas 4, 1326 (1997). 
[14] A. K. Sundaram and A. Sen, Phys. Rev. Lett. 44, 322 

(1980) . 

[15] A. Sen and A. K. Sundaram, Phys. Fluids 24, 1303 

(1981) . 

[16] A. Thyagaraja, R. D. Hazeltine, and A. Y. Aydemir, 

Phys. Fluids B 4, 2733 (1992). 
[17] A. Thyagaraja and F. A. Haas, Phys. Fluids B 5, 3252 

(1993). 

[18] W. H. Matthaeus and S. L. Lamkin, Phys. Fluids 29, 
2513 (1986). 

[19] J. F. Drake, R. G. Kleva, and M. E. Mandt, Phys. Rev. 

Lett. 73, 1251 (1994). 
[20] E.-J. Kim and P. H. Diamond, Astrophys. J. 556, 1052 

(2001). 

[21] J. Aparicio, M. G. Haines, R. J. Hastie, and J. P. Wain- 

wright, Phys. Plasmas 5, 3180 (1998). 
[22] J. D. Meiss, R. D. Hazehine, P. H. Diamond, and S. M. 

Mahajan, Phys. Fluids 25, 815 (1982). 
[23] P. K. Kaw, E. J. Valeo, and P. H. Rutherford, Phys. Rev. 

Lett. 43, 1398 (1979). 
[24] A. Lazarian and E. T. Vishniac, Astrophys. J. 517, 700 



14 



(1999). 

[25] M. Dubois and A. Samain, Nucl. Fusion 20, 1101 (1980). 
[26] M. N. Bussac and R. Pellat, Phys. Rev. Lett. 59, 2650 
(1987). 

[27] Y. Nishimura, J. D. Callen, and C. C. Hegna, Pliys. Plas- 
mas 6, 4685 (1999). 

[28] F. Porcelli, D. Boucher, and M. N. Rosenbluth, Plasma 
Phys. Control. Fusion 38, 2163 (1996). 

[29] H. R. Strauss, Phys. Fluids 19, 134 (1976). 

[30] K. Nishikawa and M. Wakatani, Plasma Physics 
(Springer, Berlin, 2000). 

[31] A. L. M. Rogister, R. Singh, and A. Kaleck, Nucl. Fusion 
29, 1175 (1989). 

[32] A. Rogister, A. Kaleck, M. Psimopoulos, and G. Hassel- 
berg, Phys. Fluids B 2, 953 (1990). 



[33] E. Westerhof, P. Smeulders, and N. L. Cardozo, Nucl. 
Fusion 29, 1056 (1989). 

[34] H. R. Koslowski, G. Fuchs, A. Kramer-Flecken, J. Rapp, 
and the TEXTOR-94 team. Plasma Phys. Control. Fu- 
sion 39, B325 (1997). 

[35] D. Biskamp and J. F. Drake, Phys. Rev. Lett. 73, 971 
(1994). 

[36] A. Thyagaraja, P. J. Knight, M. R. de Baar, G. M. D. 

Hogeweij, and E. Min, Phys. Plasmas 12, 090907 (2005). 
[37] A. J. H. Donne, J. C. van Gorkom, V. S. Udintsev, 

C. W. Domier, A. Kramer-Flecken, N. C. Luhmann, Jr., 

F. C. Schiiller, and TEXTOR team, Phys. Rev. Lett. 94, 

085001 (2005). 



